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Abstract 

We present a remarkable discretization of the classical Kepler prob- 
lem which preserves its trajectories and all integrals of motion. The 
points of any discrete orbit belong to an appropriate continuous tra- 
jectory. 
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In last years there is a growing interest in the subject of the geomet- 
ric integration of ordinary differential equations pQ. Geometric integration 
consists in numerical solution of differential equations while preserving some 
physical or mathematical properties of the system exactly (i.e., up to round- 
off error) [2]. In this Letter we propose a modification of this technique 
focusing not only on properties of the continuous system but also on meth- 
ods of generating its exact solutions. Similar approach is well known in the 
case of integrable systems of nonlinear partial differential equations [HUB]. 

In a recent paper [3] we discuss discretizations of the classical harmonic 
oscillator equation x = —x. One of these discretizations, namely 

•Kn+Y 2x n -|- X n — i . , 

4sin 2 (V2) = ~ X n K) 
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turns out to be exact. Indeed, the general solution of (JTJ 

x n = Xo cosinh) H — — sm(nh) (2) 

sm h 

satisfies x n = x(hn), where the constant h is the time step and x(t) is the 
solution of the equation x = —x with the following initial data: x(0) = xq, 
x(0) = (xi — xq cos K)j sin h. It means that the points of the discrete solution 
x n coincide with the continuous solution x(t) evaluated at the discrete time- 
lattice, and this equality is strictly satisfied for any n. 

The existence of a discretization of this kind seems to be an exceptional 
phenomenon. The list of differential equations for which exact discretizations 
were found is rather short jH] , although includes all ordinary differential equa- 
tions with constant coefficients [HI El- Actually, it is not easy even to find 
discretizations preserving the energy integral for one-dimensional Newton 
equations, however in this case the problem of finding equations admitting 
"integrable" discretizations is to some extent solved [S]. The algorithms 
preserving both the symplectic form and the energy integral are often char- 
acterized by variable time stepping [5J [TO] • 

In this Letter we consider the classical Kepler problem 

dp kr _ dr 

-dl = -T^ p = m Tf (3) 

We recall that the angular momentum L, the total energy E and the Runge- 
Lenz vector A (pointing at the perihelium), namely 

L = rxp, E = ^- k -, A = i^-k r :, (4) 
2m r m r 

are integrals of motion for (J3J). The exact solution r(t) is not known (except 
some particular cases, e.g., circular orbits) but all trajectories (orbits) can be 
exactly found. In polar coordinates all orbits are described by the formula 



p 12 L 2EL2 



1 + e cos((p — (fo) ' km ' V rnk 2 

where (fo is an arbitrary constant (the initial condition). 

We are going to solve the following problem: find a discretization of 
(0) which preserves all integrals of motion and, if possible, preserves the 
trajectories. 
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The existence of such discretization is far from being obvious. On the 
contrary, in numerical approaches to the Kepler problem it was difficult to 
preserve all the integrals of motion, not saying about the orbits ^T] . Actually, 
it is not easy to preserve even the energy integral [TT) IT2] . 

Fortunatelly, a discretization of this kind exists. What is more, this 
algorithm is surprisingly simple: 

Ap n kr n +i _, Ar n 

77 = 2 7 > Pn ■■= m ^rr , (6) 

ar~ +l r n cos o At 

where 

Ar n = r n+ \ — r n , Ap n = p n +\ — p n , At n = t n +i — t n , 

and a and S are constant parameters. The discrete time lattice t n is chosen in 
such a way that the angle A n between r n +i and r n does not depend on n. We 
denote the half of this angle by 8 n . Therefore, by assumption, 5 n = 5 = const. 
Note that A n = A = 25 and r n ■ r n+ i = r n r n+ \ cos A. As initial data for 
© we take: r , r\ and At . The time step At n can be computed using 
the condition 5 n = const, compare (|12|). In the continuum limit 5^0. 
Therefore, comparing (jSJ) and (jHJ), we see that in this limit a — > 1. 

The discrete equation © has several remarkable properties. They are 
very unusual as far as discrete analogues of continuous equations are con- 
cerned. This discretization preserves all important properties of the Kepler 
problem (JHJ), including its integral of motions and its trajectories: elliptic, hy- 
perbolic and parabolic orbits. The presence of the parameter a, not spoiling 
advantages of this discretization, can be used to obtain a better simulation 
of the continuous solutions by the discrete solutions. In this way we can fit, 
for instance, the period of the elliptical orbit. 

The integrals of motion for the discrete Kepler problem are given by: 

f 5 _ „ (Pn) 2 k x Pn*L n kR n 

L n = aR n x p n , h n = — — , A n 



2m aR n m R n 

where 

5 r n+l r n + r n r n+l D \ & \ ^ r n r n+l cos , . 

K n :— — , ti n '■— \-tin\ — 7" ■ [') 

r n + r n+l r n + r n+l 

The vector R n has a simple geometric interpretation: the end of R n lies in 
the center between the ends of f n and f n +i (in particular, R n bisects the 
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angle between r n and r n+1 ). The angular momentum can also be expressed 
directly in terms of r n and f n+ i 

-f amr n x r n+ i . 

L n Ar n = 2mar n+ ir n sin o cos o . (8) 



One can show by straightforward calculation that L ni E n and A n do not 
depend on n. Indeed, the discrete Kepler equation © is equivalent to 

?n+i + r n _i = ( 1 | __L fcAt w _i \ ^ _ 

At n At„_i yAt„ At n _! amr2r n _! cos 5 J 

Therefore, 

r f -> \ r n+l . \ „ 

L n - L n _! = amr„ x — — + — = . 

\ L\t n At n _i J 

Applying the second equation of (jS)) to L n = L n _i we get the identity 

r n+ iAt n _! = r n _iAt n . (10) 

Then, using © and (fTUj) . we compute £7 n — -E„_i = and A n — A ra _i = 0. 
All integrals of motion can be expressed by initial data: r*o, f\ and Ato- 

Now, we proceed to showing that the equation yields an explicit nu- 
merical scheme to produce f n from the initial data r , f\ and At - The crucial 
point is to find an appropriate iterative procedure for At n . The equation (fTU|) 
is not sufficient because it contains r n+1 . We have to use the constancy of 
5 n . Namely, 5 n = 5 implies (after elementary geometric considerations) 

r n+l r n-l n x r n m \ 

1 = zcoso — . (11) 

'"n+l '"n— 1 '"n 

We substitute (fTTj) into Q, eliminate r n+ i using (JTUJ), and as a result we 
obtain 

/ 2r n _i cos 5 

\ At n „! At n At n _! mr n r n _iacos5) r n 
and it yields 

n ~ r r n _ x ^(Atp) 2 ' 

2 cos d 1 H 5-5 - 

r n mrfroacoso 
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where we took into account r n r n _iAt = r 1 r At„_ 1 which follows from the 
angular momentum conservation law L„_i = L . Therefore, the complete 
explicit procedure of iterative computation of f n consists of two discrete 
equations 0, (|T2j) and initial data r , r*i, At - 

In order to obtain trajectories for the continuous Kepler problem one has 
to express the first two conservation laws (@J) in polar coordinates. The same 
approach is effective also in the discrete case. First, using (jHJ) to eliminate 
At n , we transform the expression for the kinetic energy obtaining 

(p n f _ C 2 ( r n+1 - r n \ 2 C 2 



2m m 2 a 2 (r n+ ir n ) 2 \2 sin 5 cos 8 J m 2 a 2 r n r n+ i cos 2 5 ' 

where £ = \L n \. Then, still in a direct analogy with the continuous case, we 
introduce a new variable u n : 

1 akm cos 5 . s 

Un '' = T n 7? — ' ( 3) 

and the energy conservation law can be rewritten as 

(u n+ i — u n ) 2 2m£a 2 cos 2 5 k 2 a 2 m 2 cos 2 S . s 
(2sin^) 2 + Un+lUn = T 2 + ' (M) 

where £ = E n and the right hand side of (|14p. although complicated, is 
constant (does not depend on n). Sometimes it is convenient to use another, 
equivalent, form of (JHJ): 

u n+ i — u n cos A\ 2 2 m 2 k 2 a 2 ( ' 2£C 2 \ 

' U n = "FT— 1 + — TT ■ (15) 



sinA / C 4 \ mk 2 / 

Both (|14|) and (JT5j) can be interpreted as the energy integral for harmonic 
oscillator equation. Indeed, (fTlj) implies 

fo ■ x\2 + u n = , 16 

(2 sind) 2 

which can be easily solved (compare (JTJ), ©): 

Ui -n cosA . 

M n = M cosnAH smnA , (17) 

sin A 

and, after an elementary trigonometric transformation, we get 

u n = ^ cos(nA - 6»o ) (18) 



where 



Mi— u cos A £ 2 / 2££ 2 

tan# = ^-r — , V = - , e = \l + — — . 

Mo sin A kma V ink 

From (|13|) and (|18|) we obtain general expression for the orbits in the discrete 
case: 

V 



cos 5 + e cos(n A — 9n 



(19) 



This formula is similar to (J5J). However, in order to obtain the exact dis- 
cretization of the trajectory we have to require: V = pcos5 and e = ecosS, 
which implies 

E cos 5 mk 2 sin 2 5 



£ = LVa cos 5 , £ = — — . 20 

a laL i cos o 

In the case of a periodic motion (negative E) the parameter a can be 
fixed by the requirement that the period of the motion in the discrete case 
is the same as in the continuous case. 

After this work was completed, I become aware of a sequence of papers 
by Minesaki and Nakamura on discretizations preserving integrals of motion 
O EH O El- In particular, they also succeeded to preserve all integrals 
of motion and the orbits in the case of the Kepler problem. However, the 
method used in my paper is different from their approach. It would be 
interesting to compare results generated by both discretizations. 

Acknowledgments. The work was partially supported by the KBN grant 
No. 1 P03B 017 28. The problem of finding the best discretization of the 
Kepler problem aroused in the framework of the cooperation with Boguslaw 
Ratkiewicz. 



References 

[1] E.Hairer, C.Lubich, G. Wanner: Geometric numerical integration: structure-preser- 
ving algorithms for ordinary differential equations, Springer, Berlin 2002. 

[2] R.I.McLachlan, G.R.W.Quispel: "Geometric integrators for ODEs", J. Phys. A: 
Math. Gen. 39 (2006) 5251-5285. 

[3] B.M.Herbst, M.J.Ablowitz: "Numerically induced chaos in the nonlinear Schrodinger 
equation", Phys. Rev. Lett. 62 (1989) 2065-2068. 

[4] M.J.Ablowitz, B.M.Herbst, C.M.Schober: "Discretizations, integrable systems and 
computation", J. Phys. A: Math. Gen. 34 (2001) 10671-10693. 

6 



J.L.Cicsliriski, B.Ratkiewicz: "On simulations of the classical harmonic oscillator 
equation by difference equations" , preprint physics/05071 82 (2005); Advances in Dif- 
ference Equations 2006 (2006), Article ID 40171. 

R.P.Agarwal: Difference equations and inequalities (Chapter 3), Marcel Dekker, New 
York 2000. 

J.G.Reid: Linear system fundamentals, continuous and discrete, classic amd modern, 
McGraw-Hill, New York 1983. 

Yu.B.Suris: The problem of integrable discretization: Hamiltonian approach (Chapter 
20), Birkhauser, Basel 2003. 

T.D.Lee: "Difference equations and conservation laws", J. Stat. Phys. 46 (1987) 
843-860. 

C.Kane, J.Marsden, M.Ortiz: "Symplectic-cnergy-momentum preserving variational 
integrators", J. Math. Phys. 40 (1999) 3353-3371. 

B.Leimkuhler: "Reversible adaptive regularization: perturbed Kepler motion and 
classical atomic trajectories", Phil. Trans. R. Soc. London A 357 (1999) 1101-1134. 
A.Iserles: "Insight, not just numbers", Proc. 15th IMACS World Congress, vol. II, 
cd. by A.Sydow, pp. 589-594; Wissenschaft & Technik Verlag, Berlin 1997. 
Y.Minesaki, Y.Nakamura: "A new discretization of the Kepler motion which con- 
serves the Runge-Lenz vector", Phys. Lett. A 306 (2002) 127-133. 
Y.Minesaki, Y.Nakamura: "A new conservative numerical integration algorithm for 
the three-dimensional Kepler motion based on the Kustaanheimo-Steifel regulariza- 
tion theory", Phys. Lett. A 324 (2004) 282-292. 

T.Inoue, Y.Minesaki: "A numerical integrator for the two-fixed-centers problem con- 
serving all constants of motion", J. Phys. A: Math. Gen. 39 (2006) 9437-9452. 
Y.Minesaki, Y.Nakamura: "New numerical integrator for the Stackel system con- 
serving the same number of constants of motion as the degree of freedom" , J. Phys. 
A: Math. Gen. 39 (2006) 9453-9476. 



7 



